#include "gaussLegendre.h"
#include "utils.h"

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <gmp.h>

pthread_mutex_t R;
pthread_mutex_t S;

//-----------------------------------------------
// Sequencial
//-----------------------------------------------

void gaussLegendre_sequencial (long int n) {

    long int i;
    mpf_t pi;
    mpf_t a0, b0, p0, t0;
    mpf_t a1, b1, p1, t1;

#ifdef CALCULATE_TIME
    double time = clock();
#endif

    // - - - Inicializando - - -

    // a0 = 1 |  b0 = 1/sqrt(2) | t0 = 1/4 | p0 = 1
    mpf_init_set_ui (a0, 1);
    mpf_init_set_ui (b0, 2);
    mpf_sqrt        (b0, b0);
    mpf_ui_div      (b0, 1, b0);
    mpf_init_set_d  (t0, 0.25);
    mpf_init_set_ui (p0, 1);

    // a0 = b1 = t1 = p1 = pi = 0
    mpf_init (a1);
    mpf_init (b1);
    mpf_init (t1);
    mpf_init (p1);
    mpf_init (pi);

    // - - - Caculando - - -

    for (i = 0; i < n; ++i) {

        // a1 = (a0 + b0)/2
        mpf_add    (a1, a0, b0);
        mpf_div_ui (a1, a1, 2);

        // b1 = sqrt (a0*b0)
        mpf_mul  (b1, a0, b0);
        mpf_sqrt (b1, b1);

        // t1 = t0 - p0(a0-a1)^2
        mpf_sub (t1, a0, a1);
        mpf_mul (t1, t1, t1);
        mpf_mul (t1, p0, t1);
        mpf_sub (t1, t0, t1);

        // p1 = 2*p0
        mpf_mul_ui (p1, p0, 2);

        // Atualizando os valores
        mpf_set (a0, a1);
        mpf_set (b0, b1);
        mpf_set (t0, t1);
        mpf_set (p0, p1);
    }

    // pi = ((a1+b1)^2)/4*t1
    mpf_add (pi, a1, b1);
    mpf_mul (pi, pi, pi);
    mpf_mul_ui (t1, t1, 4);
    mpf_div (pi, pi, t1);

#ifdef CALCULATE_TIME
    time = (clock() - time) / CLOCKS_PER_SEC;
#endif

    gmp_printf("%.7Ff\n", pi);

#ifdef CALCULATE_TIME
    printf("Time: %lfs\n", time);
#endif

    mpf_clear (a0);
    mpf_clear (a1);
    mpf_clear (b0);
    mpf_clear (b1);
    mpf_clear (t0);
    mpf_clear (t1);
    mpf_clear (p0);
    mpf_clear (p1);
    mpf_clear (pi);
}

//-----------------------------------------------
// Paralelo
//-----------------------------------------------

struct gaussLegendreData {

    long int n;
    mpf_t a0, b0, a1, b1;
};

void * gaussLegendre_thread (void* data) {

    // - - - Inicializando - - -

    long int i;
    struct gaussLegendreData *x = (struct gaussLegendreData*) data;

    // a0 = 1 |  b0 = 1/sqrt(2)
    mpf_init_set_ui (x->a0, 1);
    mpf_init_set_ui (x->b0, 2);
    mpf_sqrt        (x->b0, x->b0);
    mpf_ui_div      (x->b0, 1, x->b0);

    // a0 = b1 = 0
    mpf_init (x->a1);
    mpf_init (x->b1);

    // - - - Caculando - - -

    for (i = 0; i < x->n; ++i) {

        pthread_mutex_lock (&S);

        // a1 = (a0 + b0)/2
        mpf_add    (x->a1, x->a0, x->b0);
        mpf_div_ui (x->a1, x->a1, 2);

        // b1 = sqrt (a0*b0)
        mpf_mul  (x->b1, x->a0, x->b0);
        mpf_sqrt (x->b1, x->b1);

        pthread_mutex_unlock (&R);
    }

    pthread_exit(NULL);

    return NULL;
}

void gaussLegendre_paralelo (long int n, int threads) {

    // - - - Inicializando - - -
    int err;
    long int i;
    pthread_t threads_v;
    struct gaussLegendreData p_gld;

    mpf_t pi;
    mpf_t p0, t0;
    mpf_t p1, t1;

    p_gld.n = n;

    // t0 = 1/4 | p0 = 1
    mpf_init_set_d  (t0, 0.25);
    mpf_init_set_ui (p0, 1);

    // t1 = p1 = pi = 0
    mpf_init (t1);
    mpf_init (p1);
    mpf_init (pi);

    pthread_mutex_lock (&R);

    #ifdef CALCULATE_TIME
    double time = clock();
    #endif

    // - - - Calculando - - -

    // - Criando thread para calcular 'a' e 'b' -
    err = pthread_create(&threads_v, NULL, gaussLegendre_thread, (void *) &p_gld);

    if (err) {
        printf("ERRO: %d\n", err);
        exit(-1);
    }

    // - Calculando p e t na thread principal -
    for (i = 0; i < n; ++i) {

        pthread_mutex_lock (&R);

        // t1 = t0 - p0(a0-a1)^2
        mpf_sub (t1, p_gld.a0, p_gld.a1);
        mpf_mul (t1, t1, t1);
        mpf_mul (t1, p0, t1);
        mpf_sub (t1, t0, t1);

        // p1 = 2*p0
        mpf_mul_ui (p1, p0, 2);

        // Atualizando os valores
        mpf_set (p_gld.a0, p_gld.a1);
        mpf_set (p_gld.b0, p_gld.b1);
        mpf_set (t0, t1);
        mpf_set (p0, p1);

        pthread_mutex_unlock (&S);
    }

    // pi = ((a1+b1)^2)/4*t1
    mpf_add (pi, p_gld.a1, p_gld.b1);
    mpf_mul (pi, pi, pi);
    mpf_mul_ui (t1, t1, 4);
    mpf_div (pi, pi, t1);

    #ifdef CALCULATE_TIME
    time = (clock() - time) / CLOCKS_PER_SEC;
    #endif

    gmp_printf("%.7Ff\n", pi);

    #ifdef CALCULATE_TIME
    printf("Time: %lfs\n", time);
    #endif

    mpf_clear (p_gld.a0);
    mpf_clear (p_gld.a1);
    mpf_clear (p_gld.b0);
    mpf_clear (p_gld.b1);
    mpf_clear (t0);
    mpf_clear (t1);
    mpf_clear (p0);
    mpf_clear (p1);
    mpf_clear (pi);
}
